rm(list=ls())
gc()
library(lme4)
library(ggplot2)
library(ggpubr)
library(dplyr)


load("data/Data model cantons and municipalities.Rda")

data_mun_obs <- data_canton_mun$data_mun_simple

#### Rescale variables
data_canton_scale_simple <- transform(data_canton_mun$data_canton_simple, 
                               Discrimination = scale(c_summary.50..chain.1),
                               Discrimination.2.5 = scale(c_summary.2.5..chain.1),
                               Discrimination.25 = scale(c_summary.25..chain.1),
                               Discrimination.75 = scale(c_summary.75..chain.1),
                               Discrimination.97.5 = scale(c_summary.97.5..chain.1),
                               diff_log_sc = scale(diff_log),
                               log_sc = log_YesShare)

data_canton_scale_dynamic <- transform(data_canton_mun$data_canton_simple, 
                                      Discrimination = scale(c_summary.50..chain.1),
                                      Discrimination.2.5 = scale(c_summary.2.5..chain.1),
                                      Discrimination.25 = scale(c_summary.25..chain.1),
                                      Discrimination.75 = scale(c_summary.75..chain.1),
                                      Discrimination.97.5 = scale(c_summary.97.5..chain.1),
                                      diff_log_sc = scale(diff_log),
                                      log_sc = log_YesShare)


# 1. Run the mixed models 

### 1.A. Cantons ###

multi_model_slope <- lmer(diff_log ~ (Discrimination|number_mun), data_canton_scale_simple)
multi_model_slope.2.5 <- lmer(diff_log ~ (Discrimination.2.5|number_mun), data_canton_scale_simple)
multi_model_slope.25 <- lmer(diff_log ~ (Discrimination.25|number_mun), data_canton_scale_simple)
multi_model_slope.75 <- lmer(diff_log ~ (Discrimination.75|number_mun), data_canton_scale_simple)
multi_model_slope.97.5 <- lmer(diff_log ~ (Discrimination.97.5|number_mun), data_canton_scale_simple)

## Extact random slope

b <- broom.mixed::tidy(multi_model_slope, effects = "ran_vals", conf.int = TRUE)
b.2.5 <- broom.mixed::tidy(multi_model_slope.2.5, effects = "ran_vals", conf.int = TRUE)
b.25 <- broom.mixed::tidy(multi_model_slope.25, effects = "ran_vals", conf.int = TRUE)
b.75 <- broom.mixed::tidy(multi_model_slope.75, effects = "ran_vals", conf.int = TRUE)
b.97.5 <- broom.mixed::tidy(multi_model_slope.97.5, effects = "ran_vals", conf.int = TRUE)

data_canton_results <- as.data.frame(list(rand = b[b$term=="Discrimination",]$estimate,
                                          low = b[b$term=="Discrimination",]$conf.low,
                                          high = b[b$term=="Discrimination",]$conf.high,
                                          rand.2.5 = b.2.5[b.2.5$term=="Discrimination.2.5",]$estimate,
                                          low.2.5 = b.2.5[b.2.5$term=="Discrimination.2.5",]$conf.low,
                                          high.2.5 = b.2.5[b.2.5$term=="Discrimination.2.5",]$conf.high,
                                          rand.25 = b.25[b.25$term=="Discrimination.25",]$estimate,
                                          low.25 = b.25[b.25$term=="Discrimination.25",]$conf.low,
                                          high.25 = b.25[b.25$term=="Discrimination.25",]$conf.high,
                                          rand.75 = b.75[b.75$term=="Discrimination.75",]$estimate,
                                          low.75 = b.75[b.75$term=="Discrimination.75",]$conf.low,
                                          high.75 = b.75[b.75$term=="Discrimination.75",]$conf.high,
                                          rand.97.5 = b.97.5[b.97.5$term=="Discrimination.97.5",]$estimate,
                                          low.97.5 = b.97.5[b.97.5$term=="Discrimination.97.5",]$conf.low,
                                          high.97.5 = b.97.5[b.97.5$term=="Discrimination.97.5",]$conf.high,
                                          number_mun = b[b$term=="Discrimination",]$level))


###### Figure validation with credible intervals 


multi_model_time_slope <- lmer(diff_log ~ (Discrimination|number_mun:legislatur), data_canton_scale_dynamic)

b <- broom.mixed::tidy(multi_model_time_slope, effects = "ran_vals", conf.int = TRUE)

data_canton_results_time <- as.data.frame(list(rand = b[b$term=="Discrimination",]$estimate,
                                               low = b[b$term=="Discrimination",]$conf.low,
                                               high = b[b$term=="Discrimination",]$conf.high,
                                               number_mun = substr(b[b$term=="Discrimination",]$level, 1, nchar(b[b$term=="Discrimination",]$level) -3),
                                               t = substr(b[b$term == "Discrimination",]$level,nchar(b[b$term=="Discrimination",]$level) -1, nchar(b[b$term=="Discrimination",]$level))))






# Municipalities


####### Municipalities ###########


#### Rescale variables
data_mun_scale_simple <- transform(data_canton_mun$data_mun_simple, 
                            Discrimination = scale(c_summary.50..chain.1),
                            Discrimination.2.5 = scale(c_summary.2.5..chain.1),
                            Discrimination.25 = scale(c_summary.25..chain.1),
                            Discrimination.75 = scale(c_summary.75..chain.1),
                            Discrimination.97.5 = scale(c_summary.97.5..chain.1),
                            diff_log_sc = scale(diff_log),
                            log_sc = log_YesShare)

data_mun_scale_dynamic <- transform(data_canton_mun$data_mun_dynamic, 
                            Discrimination = scale(c_summary.50..chain.1),
                            diff_log_sc = scale(diff_log),
                            log_sc = log_YesShare)



multi_model_slope <- lmer(diff_log ~ (Discrimination|number_mun), data_mun_scale_simple)
multi_model_slope.2.5 <- lmer(diff_log ~ (Discrimination.2.5|number_mun), data_mun_scale_simple)
multi_model_slope.25 <- lmer(diff_log ~ (Discrimination.25|number_mun), data_mun_scale_simple)
multi_model_slope.75 <- lmer(diff_log ~ (Discrimination.75|number_mun), data_mun_scale_simple)
multi_model_slope.97.5 <- lmer(diff_log ~ (Discrimination.97.5|number_mun), data_mun_scale_simple)

b <- broom.mixed::tidy(multi_model_slope, effects = "ran_vals", conf.int = TRUE)
b.2.5 <- broom.mixed::tidy(multi_model_slope.2.5, effects = "ran_vals", conf.int = TRUE)
b.25 <- broom.mixed::tidy(multi_model_slope.25, effects = "ran_vals", conf.int = TRUE)
b.75 <- broom.mixed::tidy(multi_model_slope.75, effects = "ran_vals", conf.int = TRUE)
b.97.5 <- broom.mixed::tidy(multi_model_slope.97.5, effects = "ran_vals", conf.int = TRUE)

data_mun_results <- as.data.frame(list(rand = b[b$term=="Discrimination",]$estimate,
                                          low = b[b$term=="Discrimination",]$conf.low,
                                          high = b[b$term=="Discrimination",]$conf.high,
                                          rand.2.5 = b.2.5[b.2.5$term=="Discrimination.2.5",]$estimate,
                                          low.2.5 = b.2.5[b.2.5$term=="Discrimination.2.5",]$conf.low,
                                          high.2.5 = b.2.5[b.2.5$term=="Discrimination.2.5",]$conf.high,
                                          rand.25 = b.25[b.25$term=="Discrimination.25",]$estimate,
                                          low.25 = b.25[b.25$term=="Discrimination.25",]$conf.low,
                                          high.25 = b.25[b.25$term=="Discrimination.25",]$conf.high,
                                          rand.75 = b.75[b.75$term=="Discrimination.75",]$estimate,
                                          low.75 = b.75[b.75$term=="Discrimination.75",]$conf.low,
                                          high.75 = b.75[b.75$term=="Discrimination.75",]$conf.high,
                                          rand.97.5 = b.97.5[b.97.5$term=="Discrimination.97.5",]$estimate,
                                          low.97.5 = b.97.5[b.97.5$term=="Discrimination.97.5",]$conf.low,
                                          high.97.5 = b.97.5[b.97.5$term=="Discrimination.97.5",]$conf.high,
                                          number_mun = b[b$term=="Discrimination",]$level))


###### Figure validation with credible intervals 


multi_model_time_mun <- lmer(diff_log_sc ~ (Discrimination|number_mun:legislatur), data_mun_scale_dynamic)

b <- broom.mixed::tidy(multi_model_time_mun, effects = "ran_vals", conf.int = TRUE)

data_mun_results_time <- as.data.frame(list(rand = b[b$term=="Discrimination",]$estimate,
                                            low = b[b$term=="Discrimination",]$conf.low,
                                            high = b[b$term=="Discrimination",]$conf.high,
                                            number_mun = substr(b[b$term=="Discrimination",]$level, 1, nchar(b[b$term=="Discrimination",]$level) - 3),
                                            t = substr(b[b$term=="Discrimination",]$level, nchar(b[b$term=="Discrimination",]$level) -1, nchar(b[b$term=="Discrimination",]$level))))



data_results_canton_mun <- list(data_canton_results=data_canton_results,
                                data_canton_results_time=data_canton_results_time, 
                                data_mun_results=data_mun_results,
                                data_mun_results_time=data_mun_results_time)


save(data_results_canton_mun, file = "data/Results position municipalities and cantons.Rda")



